OCS ePS results processing + AF-BLMs (Jan 2022 update)¶
Quick attempt at AF-BLMs for OCS, based on guessed/reasonable ADMs - action starts below after setup stage). Status: basically working, may be some remaining bugs with ADMs, cos^2 and renormalisation, but trends with alignment should be correct.
For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
Additional notes below
Setup¶
In [1]:
!hostname
In [2]:
!conda env list
In [3]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
import epsproc as ep
# Plotters
import hvplot.xarray
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
plotWidth = 700
hvPlotters.setPlotters(width = plotWidth, snsStyle='whitegrid')
hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=plotWidth)) # Fix plot height! Currently not set by default.
In [4]:
# For class, above settings don't take, not sure why, something to do with namespaces/calling sequence?
# Overriding snsStyle does work however... although NOT CONSISTENTLY????
# AH, looks like ordering matters - set_style LAST (.set seems to override)
import seaborn as sns
sns.set(rc={'figure.figsize':(10,6)}) # Set figure size in inches
sns.set_context("paper")
sns.set_style("whitegrid") # Set plot style
sns.set_palette("Paired") # Set colour mapping
# Try direct fig type setting for PDF output figs
from IPython.display import set_matplotlib_formats
# set_matplotlib_formats('png', 'pdf')
set_matplotlib_formats('svg', 'pdf')
In [5]:
# xr.set_options(display_style='html')
Load data¶
In [6]:
import warnings
# warnings.filterwarnings('once') # Skip repeated numpy deprecation warnings in current build (xr15 env)
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
In [7]:
# # Scan for subdirs, based on existing routine in getFiles()
fileBase = Path('/home/paul/ePS/OCS/OCS_survey') # Data dir on Stimpy
In [8]:
# TODO: fix orb label here, currently relies on (different) fixed format
data = ePSmultiJob(fileBase, verbose = 0)
data.scanFiles()
data.jobsSummary()
In [9]:
# Fix labels for plots - currently have some mis-named files!
# Set state labels
v = ['C','B','A','X']
sDict = {n:v[k] for k,n in enumerate(range(9,13))}
for key in data.data.keys():
data.data[key]['jobNotes']['orbKey'] = key # Existing orb key, from filename
data.data[key]['jobNotes']['orbGroup'] = int(key.strip('orb')) + 1 # Corrected orb group
# data.data[key]['jobNotes']['orbLabel'] = f"HOMO - {12 - int(key.strip('orb')) - 1}"
data.data[key]['jobNotes']['orbLabel'] = sDict[data.data[key]['jobNotes']['orbGroup']]
# print(data.data[key]['jobNotes']['orbLabel'])
System properties¶
Note orbital numbering in table below:
Orbis Gamess file output orbital numbering, but energy but not grouped by degeneracy.OrbGrpis grouped numbering by degeneracy, also used by ePolyScat, and will be used for labels etc. below.
BUT - file names are 0-indexed (oops), so give OrbGrp-1 here. Sorry.
In [10]:
data.molSummary()
In [11]:
# Fix orb names! NOW FIXED ABOVE
AF-$\beta_{LM}$ (UPDATES 13/12/21)¶
13/12/21 PH
Revisiting AF cases with (hopefully) better match to experimental alignment.
- Add in cos^2(theta) metric.
- Ballpark ADMs assuming Gaussian case.
- Calculate & plot AF-BLMs for these cases.
Alignment metrics¶
Here we'll start with the molecular alignment defined in terms of axis distribution moments (ADMs), given by parameters $A^K_{Q,S}(t)$:
$ P(\Omega,t) = \sum_{K,Q,S} A^K_{Q,S}(t)D^K_{Q,S}(\Omega)$
Where $P(\Omega,t)$ is the full axis distribution probability, expanded for Euler angles $\Omega$.
This reduces to the 2D case if $S=0$:
$ P(\theta,\phi,t) = \sum_{K,Q} A^K_{Q,0}(t)D^K_{Q,0}(\Omega) = \sum_{K,Q} A^K_{Q}(t)Y_{K,Q}(\Omega)$
In the current code (v1.3.0, 16/08/21), the ADMs can be set directly, and expanded trivially for the 2D case (3D case to do). (Note that the 3D case is supported for AF-BLM calculations, but not yet for direct expansion & visualisation.)
Expectation values of $\langle cos^n(\theta, t)\rangle$ can be calculated directly from the $P(\theta,t)$ distributions:
$\langle cos^n(\theta, t)\rangle = \int_{\theta} cos^n(\theta)P(\theta,t) d\theta$
See, e.g.:
[1]
P. Hockett, “General phenomenology of ionization from aligned molecular ensembles,” New Journal of Physics, vol. 17, no. 2, p. 023069, Feb. 2015, doi: 10.1088/1367-2630/17/2/023069.
[2] Reid, K.L. (2018) ‘Accessing the molecular frame through strong-field alignment of distributions of gas phase molecules’, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 376(2115), p. 20170158. doi:10.1098/rsta.2017.0158.
Set test ADMs¶
Roughly eyeballed for a Gaussian distribution with $0<\sigma^{2}<0.2$ as per Figure 1 in Reid 2018 [2].
NOTE: not totally sure on the normalization here. Assuming othonorm sph, may be an additional sqrt(4pi) required for Y00? https://shtools.github.io/SHTOOLS/complex-spherical-harmonics.html#orthonormalized
UPDATE: setting for 4$pi$- normalised & converting to orthonorm (for numerical functions) looks good for comparison with ref. [2]. Note also additional factor of sqrt(2) in metrics at the moment - due to 0:$pi$ default domain and/or renorm(?)... should also confirm there is no additional ^2 in distribution generation.
In [12]:
# Set ADMs for increasing alignment...
tPoints = 10
# inputADMs = np.asarray([[0,0,0, *np.ones(10)], [2,0,0, *np.linspace(0,1,tPoints)], [4,0,0, *np.linspace(0,0.5,tPoints)]])
# inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)],
# [2,0,0, *np.linspace(1,2,tPoints)],
# [4,0,0, *np.linspace(0.3,2.8,tPoints)],
# [6,0,0, *np.linspace(0,3,tPoints)],
# [8,0,0, *np.linspace(0,3,tPoints)]]
# Ballpark ADMs from Fig. 1 in [2]
inputADMs = [[0,0,0, *np.ones(10) * np.sqrt(4*np.pi)],
[2,0,0, *np.linspace(0.3,2.2,tPoints)],
[4,0,0, *np.linspace(0,1.5,tPoints)],
[6,0,0, *np.linspace(0,1.5,tPoints)],
[8,0,0, *np.linspace(0,1.5,tPoints)]]
warnings.filterwarnings('ignore') # Skip repeated numpy deprecation warnings in current build (xr15 env)
ADMs = ep.setADMs(ADMs = inputADMs) # TODO WRAP TO CLASS (IF NOT ALREADY!)
# Set to separate dataset for testing
ADMs.attrs['jobLabel'] = 'Alignment data testing'
data.data['ADM'] = {'ADM': ADMs}
ep.lmPlot(ADMs, xDim='t');
Compute $P(\theta)$ distribtuions + metrics¶
In [13]:
# Axis distributions P(theta) vs. t (summed over phi)
# Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pStyle='grid', reducePhi='sum', returnFlag = True) #, facetDims=['t', 'Eke'])
Pt = data.padPlot(keys = 'ADM', dataType='ADM', Etype = 't', pType='r', pStyle='grid', reducePhi='sum', returnFlag = True) #, facetDims=['t', 'Eke']) # Real plot to check values OK.
In [14]:
# The returned data array contains the plotted values
Pt
Out[14]:
In [15]:
# Check no negative values...
print(Pt.min() > 0)
In [16]:
def computeCNs(Pt, N = None, Nrange = [2,10], Nstep = 2):
"""
Function to compute <cos^N> expectation values from P(theta) distributions.
NOTE: currently assumes 1D case. Should generalise!
"""
if N is None:
N = np.arange(Nrange[0], Nrange[1], Nstep)
# Generate cos^N basis
# cn = xr.ones_like(Pt) * (np.cos(Pt.Theta)**N) # Use multiplication to keep t coord - this fails for array-like N with broadcast errors
cn = xr.ones_like(Pt) * np.cos(Pt.Theta).expand_dims({'N':N}).T.pipe(np.power,N) # OK with transpose
Jt = np.sin(Pt.Theta) # Jacobian
norm = (Pt * Jt).sum('Theta')
cnMetric = (Pt * cn * Jt).sum('Theta')/norm
# cnMetric = Pt.dot(cn, dims='Theta') # Multiply and sum over Theta - dot version (no Jacobian or renorm)
# cnMetric = (Pt.dot(cn * Jt, dims='Theta')/Pt.dot(Jt, dims='Theta')) # Dot version with renorm
return cnMetric
In [17]:
# Compute and plot for a range of N
PtNorm = Pt #/Pt.sum('Theta') # Norm integral to 1 - now set in function
cnMetric = computeCNs(PtNorm)
# cnMetric.plot() # Basic plot
sf = np.sqrt(2) # Seem to need additional sqrt(2) normalisation? Probably SPH definition and/or integration domain and/or distribution generation params.
# sf = 1
ep.lmPlot(cnMetric*sf, xDim='t'); # With lmPlot()
In [18]:
# Seem to need additional sqrt(2) normalisation? Probably SPH definition?
# (cnMetric*np.sqrt(2)).sel(N=2).plot()
# (cnMetric*np.sqrt(2)).plot.line(x='t')
# Testing exact values
# sf = np.sqrt(2)
# sf = 2
# (cnMetric*sf).plot.line(x='t')
# (ADMs).plot.line(x='t')
(cnMetric*sf).fillna(0).hvplot.line(x='t').overlay('N')
Out[18]:
AFBLMs for aligned case¶
New plotting code to update current package versions!
In [19]:
# Compute AFBLMs for aligned case
keys = data._keysCheck(None)
keys = list(set(keys) - set(['ADM']))
data.AFBLM(keys = keys, AKQS = ADMs, selDims = {'Type':'L'})
In [20]:
# Quick stack to dict & XR dataset
dataType = 'AFBLM'
thres = 1e-2
selDims = None
Etype = 'Eke'
pDict = []
for key in keys:
subset = ep.matEleSelector(data.data[key][dataType], thres=thres, inds = selDims, dims = Etype, sq = True)
subset.name = 'BLM'
subset = ep.plotTypeSelector(subset, pType = 'r')
# pDict.append(subset.expand_dims({'Orb':[key]}))
# pDict.append(subset.expand_dims({'Orb':[data.data[key][dataType].attrs['jobLabel']]}))
pDict.append(subset.expand_dims({'Orb':[data.data[key]['jobNotes']['orbLabel']]}))
# hvDS = hv.Dataset(subset.unstack('BLM'))
# display(hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m']))
# pDict[key] = hvDS.to(hv.Curve, kdims=Etype).overlay(['l','m'])
# hv.HoloMap(pDict, kdims = ['Orb']) # Individual plots OK, but won't stack?
xrAF = xr.concat(pDict, 'Orb')
In [21]:
# hvPlotters.opts.defaults(hvPlotters.opts.Curve(height=700))
hvDS = hvPlotters.hv.Dataset(xrAF.unstack('BLM'))
hvDS.select(l=[2,4,6]).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m'])
Out[21]:
In [22]:
# With additional overlays
hvDS.select(l=[2,4,6], m=0).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['t']).layout('Orb').cols(1)
Out[22]:
In [23]:
# With additional overlays
# This seems to break legend?
hvDS.select(l=[2,4,6]).to(hvPlotters.hv.Curve, kdims=Etype).overlay(['l','m','t']).layout('Orb').cols(1)
Out[23]:
In [24]:
# Check also cross-secions
XS = np.abs(xrAF.XSrescaled) # Should be in MB
# XS = np.abs(xrAF.XSraw) # Unscaled/uncorrect XS
XS.name = 'XS/MB'
XS.hvplot.line(x='Eke').overlay('t').layout('Orb').cols(1)
Out[24]:
Versions¶
Code¶
In [25]:
# print(data.jobInfo['ePolyScat'][0])
# print('Run: ' + jobInfo['Starting'][0].split('at')[1])
In [26]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter', 'holoviews'])
Out[26]:
In [27]:
# Check current Git commit for local ePSproc version
!git -C {Path(ep.__file__).parent} branch
!git -C {Path(ep.__file__).parent} log --format="%H" -n 1
In [28]:
# Check current remote commits
!git ls-remote --heads git://github.com/phockett/ePSproc
Notes¶
07/01/22 PH
- Fixed/tested cos^2 values. (Maybe - now matching Reid 2018, but may not be correct in all cases?)
- Updated plotting with HV for BLMs (in dev stage).
- Checked XS values too - may have some -ve cases/phase issues to sort? (Cf. general AF code testing notes elsewhere.)
Versions
- TIDY version for distribution: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orb11_proc_131221-JAKE.ipynb
- Previous HOMO results: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orb11_proc_090621-JAKE.html
- Previous all-orb restuls: https://phockett.github.io/ePSdata/OCS-preliminary/OCS_orbs8-11_AFPAD_preliminary_110621.html
13/12/21 PH
STATUS: basics in place, need to trim code and fix a few things still.
Revisiting AF cases with (hopefully) better match to experimental alignment.
- Add in cos^2(theta) metric.
- Ballpark ADMs assuming Gaussian case.
TODO:
- Fix orb numering issue! (See note below.)
FIXED:
- Note 'it' and 'labels' handling changed in recent versions of ePSproc, and degenerate dim hanling now correct!
09/06/21 PH
Quick look at initial ePS results for orb 11 (HOMO), 2.5eV step size.
TODO:
- Currently AF code fails for Xarray = 0.17, issue with phaseCons setting to coord, sigh.
- AF normalisation not correct here, due to doubly-degenerate state? Summing over 'it' not working correctly either, needs a debug (or missing some settings...?).
- lmPlot() slice and cmapping to fix, in cases of missing/dropped dims, and for clims if possible.
For methods: https://epsproc.readthedocs.io/en/dev/demos/ePSproc_class_demo_161020.html
In [ ]: